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FEASIBILTY OF A BLAST WAVE 
ATTENUATION STRUCTURE 



by Dale Richard Hartmann 

Chairman of the Supervisory Committee: Professor Ashley F. Emery 
Department of Mechanical Engineering 

This thesis begins with an overview of bombings in the United States, followed by the 
introduction of the Rankine-Hugoniot equations for blast wave pressure. The subsequent 
chapters develop the one-dimensional and two-dimensional Euler equations. These 
equations are the solved using the MacCormack finite difference algorithm. The basis of 
the investigation then begins by placing pole, shear plate and wedge obstacles in the path 
of the blast wave. The results of these simulations are interpreted and conclusions 
presented. Finally a synopsis of the existing results and cost analysis for structure 
hardening are presented. 
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PREFACE 



The inspiration for this thesis was developed over the course of two years. While at an 
assignment with the Defense Nuclear Agency, 1 was privileged to participate in several 
tests involving car bombs. 1 became very interested with the effects of the blast waves on 
different types of structures. After familiarizing myself with the current methods of 
protecting an existing structure from blast waves, I thought there must be a better way. 
This thesis is my investigation regarding whether there is a better way to protect existing 
structures. 
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INTRODUCTION 



BOMBS: THE UNKNOWN MENACE 



The threat of terrorist bombings is present every year. Although the United States has 
been relatively immune from terrorist bombings such as those in Ireland or Israel, 2577 
bombings occurred in the United States in 1995 alone. Bombings involving improvised 
explosive devices such as pipe bombs numbered 1,562 in 1995. The majority of the 
bombings in 1995 were small in size but produced $105 million damage and 937 
casualties (193 killed and 744 injured). The worst bombing in the history of the United 
States, the Muirah building in Oklahoma City, accounted for $100 million damage and 
786 casualties (168 killed and 518 injured). 

A cursory review of the above facts and events leads to the question: How can buildings 
that are possible terrorist bombing targets be protected from the effects of bomb blasts? 
That question is the basis of this thesis. 



CHAPTER 1: THE PROBLEM 



TYPES OF EXPLOSIVES 

Explosives classified as either: High explosives or Low explosives. High explosives 
possess certain characteristics that are vastly different than Low explosives. High 
explosives detonate, which means that when they are initiated a shock or blast wave will 
form. They also will burst or shatter materials near them, are capable of penetrating 
materials, and have the capability of lifting or moving objects. 

Low explosives, in contrast, do not detonate but bum very rapidly. Since Low explosives 
do not detonate, the pressure rise that is produced is usually smaller in amplitude but 
longer in duration than that of a High explosive. This combination tends not to produce 
an impulse type of blast wave but a slightly ‘softer’ shock to materials nearby. 

BLAST WAVE CHARACTERISTICS 

A blast wave generated as the result of the initiation of a contained High explosive is 
created in the following sequence of events. First hot gases with temperatures of the order 
3000 degrees C are generated. In concert with this temperature rise the pressure of these 
gases is of the order 100 to 300 kilobars. It is this second characteristic of high explosive 
initiation that is of concern here. 

This hot high-pressure gas then expands into the surrounding atmosphere. As expansion 
occurs, the air surrounding the expanding gas is forced outward. Since air is a 
compressible medium, a layer of air adjacent to the outer edge of the expanding gas is 
compressed. It is this layer of compressed air, which is called a blast wave. As the hot 
gas expands and cools, its pressure falls. The pressure of the layer of compressed air also 
falls as it is pushed farther outward from the point of detonation. As the gas continues to 
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cool and expand, at some point in time and space the pressure will fall below atmospheric 
pressure due to the momentum of the layer of compressed air. This slight negative 
pressure is called overexpansion and results in a negative phase of the blast wave 
whereby the outward flow of gases and air is reversed. Eventually after some number of 
oscillations, equilibrium is reached and the motion of the air stops. 

The speed at which the blast wave travels is different depending on the medium in which 
it is traveling. For a bomb resting on the ground, two waves will exist: one traveling 
through the ground and another traveling through the air. The latter is what this thesis is 
concerned with. 

The magnitude of the impulse of an explosive blast wave is dependent upon many 
factors. Some of these factors are: 

Is the explosive cased? 

What type of explosive is used? 

The quantity of explosive used? 

What weather conditions are present at the time of detonation? 

What is the terrain? 

Are any structures nearby? 

PROTECTING EXISTING STRUCTURES 
THE TRADITIONAL APPROACH 

The traditional approach to protecting an existing structure from explosive blast damage 
is to harden it, that is to increase the structural strength of single or multiple components. 
Hardening of a building consists of several methods. 
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The first and most rudimentary method is to harden the glazing, if any is present, on the 
exterior of the structure. Even a relatively small bomb can cause large amounts of injury 
and death without structurally damaging a building. This carnage is achieved by 
shattering normal window glass and propelling the glass fragments at high speed 
throughout the interior of rooms along the outside perimeter of the building, anyone 
present in the rooms is impaled with large shards of glass. The use of tempered or 
tempered safety glass can minimize this effect. However, the use of blast curtains (heavy 
polymer curtains weighted along the lower edge) or the installation of Mylar film with 
reinforced window frames can eliminate this hazard. The problem is the curtains have to 
be kept closed at all times to provide protection, thus rendering the window useless as an 
architectural entity. In addition the Mylar is relatively expensive. 

The most complex, and best, method is to harden the entire structure against the effects of 
explosive blast waves. This generally consists of a major reconstruction of the interior 
structural elements of a building or encasing the building inside of another more robust 
exterior structure capable of withstanding the blast wave. This method, although highly 
effective, is also difficult to design, disruptive to the occupants, time consuming to 
construct and very expensive. 

THE NEW APPROACH 

The approach that this thesis investigates is the use of a blast wave attenuation structure. 
The form of the blast wave attenuation structure would be such that as the blast wave 
travels through it the energy of the wave is dissipated. Ideally, the energy of the wave 
would be lowered to such a degree as to produce and oveipressure of no more than one 
half atmosphere. This lowered pressure, although still high enough to cause some 
damage, is safe for most commercial types of construction. 

The forms of the attenuation stmcture to be investigated include shear plates, a field of 
poles, and of wedges. The success of a type of stmcture will depend upon its ability to 
attenuate the pressure in less than 15 meters, constructability, low cost, and aesthetics. 



CHAPTER 2: THE QUESTION 



MODELING BLAST WAVES 
BLAST WAVE PARAMETERS 

The modeling of a blast wave must address two issues. The first is the modeling of the 
static aspects of the explosion. These static aspects include pressure, density, 
temperature and shock wave velocity. The second, and more difficult, is the modeling of 
the dynamics of the blast wave motion and interactions with objects. These dynamic 
aspects are the same quantities as the static aspects but are concerned with how these 
quantities change with time and position. 

STATIC ASPECTS 

It is important to model the static aspects of an explosion since this is what provides the 
driving force for the dynamic solutions and simulations. The first static aspect to be 
modeled is the pressure generated by the detonation of the explosive. The pressure 
generated is dependent upon several factors: 

Type of explosive 

Distance from explosive 

Position of explosive relative to the ground 

The type of explosive is important since each different type of explosive contains a 
different amount of energy per unit mass. These differences are summarized in the table 
below. 
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Table 1 Explosive Equivalents 



Explosive 


Mass Specitic Energy 
Qx( kJ/kg) 


TNT Equivalent 
(Qx/Qtnt) 


RI)X (Cyclonite) 


5360 


1.185 


HMX 


5680 


1.256 


Nitroglycerin (Liquid) 


6700 


1.481 


TNT 


4520 


1.000 


Blasting Gelatin 


4520 


1.000 


Nitroglycerin dynamite 


2710 


0.600 


Semtex 


5660 


1.250 


Compound B 


5190 


1.148 



As can be seen from the table above, the method of normalizing the energies of different 
explosives relative to that contained in TNT has been developed and is universally 
accepted. 

The relationship between the range and the TNT equivalent charge mass is known as the 
scaled distance and is given by the following equation: 



Z = 



R 

W m 



( 1 ) 



Where W, the universal symbol for TNT equivalents, denotes the mass of explosive and 
R denotes the range from the detonation of the explosive to the point of interest. Rankine 
and Hugoniot produced the first analytic solutions for blast wave front parameters in 
1870 for normal shocks in ideal gases. These solutions were later expanded by Brode to 
determine the peak static overpressure in the near and medium fields of the blast wave. 
The near field is the region where the overpressure is higher than 10 bar ( 1 0 6 Pa). The 
medium field is the region where the overpressure is less than 10 bar (10 6 Pa) and higher 
than 0. 1 bar (10 4 Pa). The equation for the oveipressure in the near field is given by 



Ps = 



^■+l)(l00000) (Pa) 



( 2 ) 



The equation for the oveipressure in the medium field is given by 
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0.975 1.455 5.85 

— — + — ^ + — T- 



V z z 2 z 3 



-0.019 (100000) (Pa) 



(3) 



The absolute pressure of the detonation is given by 




(4) 



The equations developed by Rankine, Hugoniot and Brode, shown above, dealt with 
shock waves in free air. These shock waves were allowed to propagate in all directions, 
or, more simply, in a spherical manner. Since I am concerned with terrorist bombings, 
the majority of which involve explosives placed on or near the ground, a correction factor 
is needed on pabsoiute to account for the reflection from the ground. If the ground were an 
ideal surface for reflection this correction factor would by 2. However, the ground will 
respond in two ways that are not ideal. The ground underneath the explosion will 
compress. The area around the explosion will undergo a fracturing process, which will 
result in particles being ejected into the air. There is no analytical method of quantifying 
these responses, but empirical evidence suggests a correction factor given by 



This correction factor is applied to the pressure calculations for an explosion to provide 
the input pressure for all analytical and finite difference calculations. The Mach number 
of the detonation shock wave is given by 



Pa = IS P, 



absolute 



(5) 



M 



x 



P ambient JJ 




( 6 ) 



The units in the equations below are correct for p a expressed in Pa. 



The density of the gas behind the detonation shockwave is given by 
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P 2 = 



P ambient 



(Kg/nv ) 



y + V 



M \ 



The velocity of the detonation shockwave is given by 



r 



u 2 =ca 



y + 1 



m 

M t 



(m/s) 



Where Ci is the ambient speed of sound 



1/2 



c, =1 y— -- - ] (m/s) 

k P ambient ' 



The temperature of the gas behind the detonation shockwave is given by 



r : =-^-(K) 
P 2^ 



Where 



With 



R = 



_A_ 

MW air 



(J/kgK) 



R u = 8314.51 (J/kmole K) 
MW air = 29 (Kg/kmole) 



(7) 



( 8 ) 



(9) 



( 10 ) 



( 11 ) 



( 12 ) 

( 13 ) 



The speed of sound behind the detonation shockwave is given by 
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c , = 



1/2 



T 

ambient * 



(m/s) 



(14) 



DYNAMIC ASPECTS 



THE EULER EQUATIONS 



Air has a small thermal conductivity (k) and also a small viscosity (o) therefore the terms 
in the Navier-Stokes equations that depend on k and o can be neglected, which leads to 
equations in which the convective terms dominate and the fluid (in this case air) is treated 
as inviscid. The inviscid 2-D Navier-Stokes equations are called the Euler equations and 
are the equations dealt with in this paper. These equations are much simpler in form and 
also programming effort. 

The Euler equations are written as 



where 



dQ cF PG 
— + — + — 



dt 3c cfy 



0 



(15) 



p ' 




pu 




pv 


pu 




pu 2 + p 




puv 




F = 




G = 


2 


pv 




puv 




pv- +p 


_E 




u[E + p) 




v(E + p) 



( 16 ) 



and 
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E = p 



e + — 
V 2 



(,r + v=) 



(17) 



The Euler equations represent a system of equations that conserve mass, momentum and 
energy. These equations must be solved as a system to capture nonlinearities such as 
shocks. 

The Euler equations also need an equation of state, which is the ideal gas law. 

p = pRT = ( r -\){E-'-f{u'+v^ (18) 



CHAPTER 3: ONE DIMENSIONAL EQUATIONS 



THE 1-D EULER EQUATIONS 

The first step in modeling the blast wave is to construct a one-dimensional model with 
time and distance as the parameters. The 1-D Euler equations are written as 



where 



and 



3 3c 





p " 




pit 


Q = 


pu 


F = 


pit 2 + p 




E 




u(E + p ) 



E = A e+ 2 



The equation of state is 



p = p RT = {y- 1 )[e - 



(19) 



( 20 ) 



( 21 ) 



( 22 ) 



This system of equations yields the following three equations 
Conservation of mass 



dp | c{pu) _ Q 
3 3c 



( 23 ) 
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Conservation of momentum 



c{pu) (\pu 2 + p) 
3 3c 



( 24 ) 



Conservation of energy 



35 3i(E + p) _ 
3 3c 



( 25 ) 



FINITE DIFFERENCING 

To evaluate this system of equations the method of finite differencing is used. Finite 
differencing is a powerful numerical technique that can be used to solve partial 
differential equations. Finite differencing algorithms come in two varieties: explicit and 
implicit. An explicit algorithm solves each equation in turn and then applies those results 
to the initial values for the next time step. An implicit algorithm solves the system of 
equations simultaneously via matrix manipulations. I chose the explicit form of finite 
differencing instead of the implicit form for the following reasons: 

Ease of coding 

Lower computational requirements 
Ability to run on desktop computing systems 

The first step in constructing a finite differencing program is to choose an algorithm. 
Many different algorithms exist but to be useful for the problem at hand the chosen 
algorithm must be able to “capture” the shock discontinuities well and also handle 
surface interactions. Several algorithms meet these requirements: 

Forward Time Centered Space (FTCS) 

MacCormack 



Harten-Ye 
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1 st Order Roe 

1 initially investigated the use of the FTCS algorithm, however 1 was not pleased with the 
shock “capturing” abilities of this algorithm. This algorithm smears or spreads the shock 
discontinuity over 4 to 5 Ax intervals, it also has dispersion errors which appear as 
oscillations at the shock discontinuity. The following figure demonstrates these points. 






Figure 1 FTCS Output 



The next algorithm I investigated was MacCormack’s. This algorithm is a two-step 
method that is known for “capturing” shock discontinuities very well. In contrast to the 
FTCS algorithm, MacCormack’s does not have the dispersion errors at the shock 
discontinuity. The following figure demonstrates these points. 
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Dislanca (m) 



Figure 2 MacCormack Output 



The Harten-Yee and 1st Order Roe methods were abandoned early on due to their 
complexity and coding difficulties and because they have been shown not to be 
substantially superior to MacCormack’s method. 

Due to the simplicity and shock capturing abilities of the MacCormack algorithm, I have 
chosen to use it as the basis for the remainder of this thesis. 

The method of the MacCormack algorithm is contained in the two steps below: 

Predictor step: Qj = Q” - AtA x F” (26) 

Corrector step: Q. n+] = +£>/)- yV, F. (27) 



where 
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A, = v - f; (28) 

v *=*}-^ (29) 

For the one dimensional case the MacCormack algorithm equations become 
Conservation of mass 



dp_ | d(pu) = Q 

Ft 3c 

Predictor: J, = pf - ((/tit ) .^ - (/*),") 

Corrector: p"" = P ‘ P ‘ — ^(/w/ - pu,_," J 
Conservation of momentum 



(30) 



(31) 



(32) 






(33) 



Predictor: (pw) . = {pu)" 



At 


'(Hi’) 


~Ax 






(34) 



Corrector: (/?h )." +1 




( 35 ) 



Conservation of energy 
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cE di{E + p) 



( 36 ) 



dt dx 



Predictor: Ej = E " 




( 37 ) 



Corrector: Ej ' ' 1 




The equations for the conservation of momentum and energy above contain a pressure 
term that is evaluated using the equation of state. The equation of state is the ideal gas 
law and is given by 



In the above form the equation state can be evaluated using the results of the 
MacCormack algorithm. 

The final concern for any finite differencing algorithm is stability. The use of the 
Courant-Friedrichs-Lewy or CFL number is the most widely used method of controlling 
stability. The CFL number arises out of the results from a Von-Neumann stability 
analysis. The theoiy of Von-Neumann utilizes a Fourier transform to transform the finite 
difference solution into a wave space solution. The amplitudes of the waves in this 
resulting solution will grow or decay based upon the particular finite difference algorithm 
chosen and the value of the CFL number. The waves in the Von-Neumann solution that 
grow are unstable, those that decay are stable. The results of the Von-Neumann analysis 
require that the CFL number be less than or equal to one. The physical implication of the 
CFL number being less than or equal to one is that the sum of the amplitudes of the 
waves that decay is larger than the sum of the amplitude of the waves that gl ow. Thus if 




( 39 ) 
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the Courant number remains below the value of 1 the decaying waves will dominate and 
the algorithm will remain stable. 



CFL = 



Ay 

(|«| + c)A/ 



( 40 ) 



where c is the speed of sound for the given pressure and density and CFL is the Courant- 
Friedrichs-Lewy number. To determine the time step used in the finite differencing code 
the CFL equation above is manipulated such that the following equation for a one- 
dimensional algorithm time step results. 



A t = CFL 



Ay 

(M + c ) 



(41) 



VALIDATING THE ALGORITHM 

At this point, after developing the equations for the one-dimensional case and 
constructing the finite difference algorithm, validation of the method is appropriate. To 
validate the code a model was constructed with a right traveling shock that reflects from 
an infinite wall at the right most-boundary. I selected this validation method due to its 
ease of coding and readily available analytic solution. 



P 2 

c 2 



► 



u 



2 




Pi 

► 

C 1 



Figure 3 Shocktube with right traveling wave 



The actual validation was done by allowing the finite difference solution to run long 
enough following the interaction with the right most wall and then comparing the 



18 



pressure and density values from the model to the analytic values calculated from the 
analytic solution (Shapiro) 

The constant initial conditions in the program code are 



Table 2 Fixed initial conditions 



Temperature 


Gamma Density 


Pressure 


Speed of sound 


298 K 


1.4 1.614 kg/m J 


101325 Pa 


296 m/s 



The initial conditions input to the finite difference model were 



Table 3 Input initial conditions 



Courant 

number 


Calculation 

distance 


Number of 
position 
steps 


Number of 
time steps 


Distance 
from charge 


Explosive 


Explosive 

mass 


0.5 


4 m 


101 


250 


10 m 


TNT 


200 kg 



The program calculated values behind the shock prior to the wall interface were 
determined using the Rankine-Hugoniot and Brode equations, equations 2, 3, 5, 6, 7, 8, 9, 
10 and 14. The values output from the program are tabulated below. 



19 



Table 4 Calculated fl ee stream shock values 



Static pressure of blast 


5.7417x10° Pa 


Mach number of shock 


2.2504 


Density 


4.8729 kg/m J 


Shock velocity 


446. 1 8 m/s 


Temperature 


410.97 K 


Speed of sound 


348.97 m/s 



The calculated values from the program behind the shock after the wall interface were 



Table 5 Calculated reflection shock values 



Static pressure of blast 


2.246x10° Pa 


Density 


■ 2.179 kg/m°S 



To obtain the analytic values an iterative method was used. Following the reflection 
from the right most wall the shock is now traveling to the left. The analytic solution is to 
hold the shock stationary and iterate for the value of W. 
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u 2 



c 



2 



P 2 



W 



► 



w + u 




u, 



Ps 



Figure 4 Shocktube with stationary wave 



The equations required for the iterative solutions are 



AA W+U 2 

M = - 



£l =1+ J1L( m 2 

Pi r + 1 1 ' 




p± 

pi 



i 



l- 




r + 1 v 




2 f 



t / 3 — u 2 — c 2 



,r + 1 






M -■ 



\ 



MJJ 



(42) 



(43) 



(44) 



(45) 



Using the equations above, equations 41 through 44, with the initial values listed in the 
table below, the results of the analytic solution will be those listed in the table below. 
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Table 6 Analytic results for shocktube 



Initial values 


1‘2 


C2 


P2 


p 2 y 




446.18 m/s 


348. 1525 m/s 


574.1 7 kPa 


4.8729 kg/m 3 1.4 


Calculated values 


W 




Ps 


Ps 




180.5 m/s 




2.0747 MPa 


11.4963 kg/m 3 



Comparing these analytic values with the finite difference calculated values yields an 
error of +6% for the density and +8% for the pressure. Which means that the finite 
difference program will yield a slightly conservative (higher than true values) result. 
This is acceptable for my puiposes. 



CHAPTER 4: TWO-DIMENSIONAL EQUATIONS 



THE 2-D EULER EQUATIONS 

The addition of a second dimension to the 1-D Euler equations is fairly straightforward. 
The second dimension requires one additional equation to the matrix and two additional 
terms. The 2-D Euler equations are written 



where 



and 



dQ cF cG n 

— + + = 0 

dt 3c dy 



p " 




pu 




pv 


pu 




pu 2 + p 




puv 


F = 


G = 


2 


pv 




puv 




pv +p 


_E _ 




u(E + p ) 




v(E +p) 



E = f\e + -(ir+v'-) 



(46) 



(47) 



(48) 



The equation of state is 



p = pRT = (r-\)[^E-^f{ir +v 2 ^ 



( 49 ) 



This system of equations yields the following four equations 
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Conservation of mass 



dp_ | cjpu) | c(pv) _ 
3 3c dy 



( 50 ) 



Conservation of x momentum 



3{pu) | d(pic + p) | ojpuv) Q 
3 3c dy 



(51) 



Conservation of y momentum 



3{pu) | 3{puv) | + p) _ Q 

3 3c dy 



(52) 



Conservation of energy 



cE 3i{E + p) 3o{E + p ) 
3 3c dy 



(53) 



The method of the MacCormack algorithm for two dimensions is similar to that for one 
dimension. The two dimensional algorithm is contained in the two steps below 



Predictor step: 


e,-.» = e,/ - a/(a,v + a,<v) 


(54) 


Corrector step: Q J k " +x ~ 


+v,c w ) 


(55) 



where 



A = F " - F " 

a X r j+\,k r j,k 



(56) 



V = F - F 

y X r j,k r j-\,k 



( 57 ) 
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and 



= G jM ; - Gj/ 



Vy = G j,k - G j,k - 1 



(58) 



(59) 



For the two dimensional case the MacCormack algorithm equations become 
Conservation of mass 



Predictor: p jk 



Corrector: p ]<k " + ' 



dp_ | c\pii) | c{pv) = Q 
& 3c dy 



= Pj,k - ^ 






„\ 



Ay 



Ay 
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Conservation of x momentum 
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Conservation of y momentum 
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where 
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Conservation of energy 
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The pressure terms in the conservation of energy equations above are evaluated by using 
the equation of state. The pressure terms are distinct for each direction. The x direction 
pressure term is given by 
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Similarly the y direction pressure term is given by 
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The time step in two dimensions is a more complicated expression given by 
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where c is the speed of sound for the given density and pressure and CFL is the Courant- 
Friedrichs-Lewy number. 

VALIDATING THE ALGORITHM 

Applying reasoning similar to the one-dimensional case, an input pressure was applied to 
the x edge of the computational domain. The effect of this input pressure was the same 
as the calculations preformed by the one-dimensional code. Inputting the same 
parameters into both the 2-D and 1-D programs allowed comparison of the output 
matrices, which were identical. Therefore the 2-D code is validated since the 1-D code 
has been proved correct and acceptable. Following validation of the 2-D code for a linear 
pressure wave traveling in the y direction, the same procedure was applied for a wave 
traveling in the x direction. The same results were obtained. 

The results of this validation are shown in the tables below 



Table 7 1-D program calculated values 



Density 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


Momentum 


2174.2 


974.5 


79.3 


0 


0 


Energy 


1.9395e6 


1.0087e6 


0.3 144e6 


0.2533e6 


0.2533e6 



The results for a 2-D wave traveling in the x direction are given in the tables below 
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Table 8 2-D x-traveling wave program calculated 
values 



Density 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


4.8729 


3.0026 


1.7713 


1.6140 


1.6140 


Momentum 


2174.2 


974.5 


79.3 


0 


0 


2174.2 


974.5 


79.3 


0 


0 


2174.2 


974.5 


79.3 


0 


0 


2174.2 


974.5 


79.3 


0 


0 


2174.2 


974.5 


79.3 


0 


0 


Energy 


1 9395e6 


1.0067e6 


0.3144e6 


0.2533e6 


0.2533e6 


1.9395e6 


1.0067e6 


0.3144e6 


0.2533e6 


0.2533e6 


1.9395e6 


1.0067e6 


0.3144e6 


0.2533e6 


0.2533e6 


1 9395e6 


1.0067e6 


0.3144e6 


0.2533e6 


0.2533e6 


1.9395e6 


1.0067e6 


0.3 144e6 


0.2533e6 


0.2533e6 



The results for a y traveling wave are in the tables below 
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Table 9 2-D y-traveling wave program calculated 
values 



Density 


4.8729 


4.8729 


4.8729 


4.8729 


4.8729 


3.0026 


3.0026 


3.0026 


3.0026 


3.0026 


1.7713 


1.7713 


1.7713 


1.7713 


1.7713 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


1.6140 


Momentum 


2174.2 


2174.2 


2174.2 


2174.2 


2174.2 


974.5 


974.5 


974.5 


974.5 


974.5 


79.3 


79.3 


79.3 


79.3 


79.3 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


Energy 


1.9395e6 


1.9395e6 


1.9395e6 


1.9395e6 


1.9395e6 


1.0067e6 


1.0067e6 


1.0067e6 


1.0067e6 


1 0067e6 


0.3144e6 


0.3 144e6 


0.3144e6 


0.3144e6 


0.3144e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 


0.2533e6 



The results above were achieved with the following input conditions: 



Table 10 2-D validation initial conditions 



Time step 


x-grid 


Range to 

Charge 


Explosive Type 


Explosive 

Mass 


5e-4 s 


5 


10m 


TNT 


200 kg 



CHAPTER 5: TESTING THE NEW APPROACH 



The validation of the 2-D code, presented in the previous chapter, was performed with 
linear wavefronts traveling in only one direction. While this situation might occur with a 
very large explosion at a great distance, in general this is an unrealistic case. Therefore, 
the use of symmetry was applied. The final step in the code development process was to 
place the explosion in one comer of the computational realm and allow the wave to 
expand in true 2-D fashion. Placement of the explosive in this corner allows for two 
planes of symmetry and reduces the number of calculations per time step by a factor of 
two. This reduction of calculations is also beneficial since it ignores the blast wave that 
travels away from the obstacles. 

The initial parameters for an explosion were placed into one comer of the computational 
realm and allowed to expand freely with no obstacles present. The results obtained were 
consistent with the expected outcome. This step was necessary to provide a test of the 2- 
D code in the true two dimensional environment. 

OBSTACLE BOUNDARY CONDITIONS 

To obtain meaningful results with obstacles in the path of the wave, a thorough 
understanding of the correct boundary conditions for the obstacles is needed. The 
minimum size of an obstacle is limited to a grid of 3 by 3 computational points. This 
minimum size is necessaiy to allow for the specification of conditions inside the obstacle. 
There is no limit on the maximum obstacle size, although it obviously cannot exceed the 
size of the computational realm. The figure below illustrates a minimum size obstacle. 
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Figure 5 Obstacle boundary conditions 



As can be seen in the figure above a 3 x 3 grid is the smallest obstacle possible to provide 
for the specification of quantities inside an obstacle. The figure above also shows what 
quantities exist at the points on the boundaries. The quantities that exist at the interior 
point are density and energy. 

For inviscid flow on rigid surfaces four distinct boundaiy conditions exist. The first is 
that the component of the x-momentum normal to the surface is zero. Similarly, the 
component of the y-momentum normal to the surface is also zero. The third condition is 
that the pressure gradient is zero. 
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The last boundary condition is that velocity component normal to the obstacle surface is 
zero. 

Another assumption needed to evaluate all parameters is constant energy inside the 
obstacle, E. The equation for energy is: 



Since at this point the velocities are known, the energy and the pressure can be 
calculated. This equation along with the velocities allows for calculating p on the 
surface. 

The pressure at the surface of the body obtained by solving the ideal gas law equation, 
which is: 




( 79 ) 




( 80 ) 
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RESULTS OF VARIOUS ATTENUATION STRUCTURES 

The following sections detail the results of running the simulation with poles, shear plates 
and wedges for obstacles. 

Before reviewing the results of the simulations it is helpful to examine the flow field 
results without any obstacles present. All of the results to be presented are at 30 time 
steps. The program does utilize variable time stepping so the results are not at the same 
total time. These results are shown in the figure below. 




Absolute Pressure 



2 4 6 6 10 12 14 16 

X grid points 




Figure 6 2-D flow field results without obstacles 

The figure above illustrates two points unique to numerical solutions. The first point is 
the gradient at the shockffont. The gradient of the shock is shown as the lines spaced 
closely together. These lines follow the curved path of the shockwave and travel at the 
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speed of the shock. The other point that is well illustrated, are the numerical oscillations 
inherent in this type of solution. The regions outlined behind the shockfront are, in fact, 
numerical oscillations of the solution and do not exist in the real explosion 



POLES 

The results for the field of poles were disappointing since it was the conceptual basis for 
this thesis. The pole obstacles, in the figure below, were set at 3 by 3 grid points with 
infinite height. The actual physical size of the obstacles varied with the total number of 
grid points and the calculation distance. For example a 10 by 10 grid, with a 4 meter 
calculation distance generates a pole that is 1.2 maters on a side. It is very important to 
ensure that the obstacle size is reasonable. In the above example to make the pole 
dimensions smaller a greater total number of grid points would be needed or a smaller 
calculation distance. 

Following the completion of many simulations the effect of poles on reducing the 
pressure of a blast wave it was found that the effect is minimal. The blast wave does 
compress the air on the blast side of the pole as 1 expected, but the shadow of the pole 
does not extend far enough in space to actually cause a decrease in the pressure of the 
wave after it has passed through the entire field. The cause of this phenomenon is the 
low viscosity of air. The calculations that are used for the simulation are inviscid, which 
allows the air to flow easily around a small object such as a pole. For the effect of an 
obstacle to be maintained farther downrange, the obstacle must be of a size large enough 
to produce a significant wake or a field of small poles set in a grid or staggered grid 
pattern that are close together. The field of poles would need to be set on the order of 
their diameter to be effective. Due to simulations running on a desktop computer system 
I could not generate a field of poles of sufficient density to demonstrate this directly. I 
am extrapolating a solution based upon the limited observations and conclusions I made 
following the simulations. The following figure illustrates this result. 
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Figure 7 Flow field results for pole obstacles 

This figure shows that the first pole, located at grid point (3,3) has a very large pressure 
gradient on its face. This gradient is produce by the compression of the air between the 
immovable pole and the moving air behind. This compressed air contains a great deal of 
energy that remains stagnant upon the face of this pole. 

The figure also shows that following passage of the shock through the pole field the 
gradient has been spread over a larger area, since the contour lines are farther apart, but 
the magnitude of the gradient has not changed. Therefore, the shock is not as steep as it 
was originally but is of essentially the same magnitude. Thus, the field of poles has 
‘softened’ the shock but has not attenuated the peak pressure. The interpretation of this 
result is that this spacing of poles, which is an example of many simulation runs, does not 
achieve the desired outcome. Extending this interpretation further leads to the conclusion 
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that no economical pole field would attenuate a shock to a level low enough to protect a 
structure. 

SHEAR PLATES 

A shear plate is a thin long plate with the long axis ideally placed parallel to the direction 
of flow. Due to the uncertainty in the location of an explosive device, I placed the shear 
plates such that the length axis is perpendicular to the face of the protected building. The 
shear plates, in the figure below, are 3 grid points wide by 5 grid points long. I 
performed numerous simulations with other sizes and orientations but this combination of 
width, length and placement best demonstrates the effect of shear plates on the shock. 

The shear plates performed slightly better than the pole field largely due to their size and 
orientation. The first plate redirects the majority of the flow field in such a way that it no 
longer raises the pressure on the lower surfaces of subsequent obstacles. The pressure on 
the upper side of each shear plate is significantly less than the pressure on the lower side. 
This is the result of the lower plates acting as walls and shielding the upper plates from 
the blast pressure. The plates also redirect the flow by straightening it in the length 
direction of the plates. After the flow has passed the end of each plate it begins to 
diffract. The placement of the shear plates is critical to achieve this effect. 



The results of the simulation can be seen in the following figure. 
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Figure 8 Flow field results for shear plate 
obstacles 

The figure above demonstrates the effect of shear plates on the flow. The lower plate 
redirects the majority of the flow along the lower edge of the computational realm and 
removes a great deal of the energy from the flow that impacts the middle and top plates. 
The flow then diffracts around the lower plate after it has passed. The gradient is again 
‘softened’, but its magnitude remains essentially the same as the original unimpeded 
shock. The large gradient in the lower left comer is produced by the shock impacting 
upon the face of the lower shear plate, which is immovable. This impact compresses the 
air and produces a large pressure peak in a very small area. The middle and top shear 
plates have similar gradients on their faces but the magnitude is much smaller. 

To provide significant protection for a building the location of the explosive would either 
need to be known or surmised before construction of the plates. The plates could then be 
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placed in such a way as to direct the flow away from the protected structure. 
Unfortunately, foresight is rarely this accurate and the construction of many plates (or 
walls for larger buildings) will be neither aesthetically pleasing nor particularly 
inexpensive. 

WEDGES 

A wedge is to two thin plates placed at an angle to each other and joined at the apex. The 
wedges used in this simulation consist of two plates 3 grid points by 5 grid points placed 
perpendicular to each other. The wedges were then placed into the computational realm 
such that the lower face of the wedge is diagonally across the path of the shock. The 
wedge in this position acts somewhat as a shield or wall in the path of the shock. 

The results for the wedge obstacles were similar to the shear plates. A large pressure 
gradient builds on the explosion side of the lower face and a shadow, or wake forms 
behind the wedge in the hollow interior. As the flow progresses past the end of the 
wedge diffraction begins and the shadow disappears. The main difference between the 
shear plates and the wedge obstacle is that, due to its size and geometry, the wedge does a 
better job of lowering the pressure in its wake. The wake, however, disappears rapidly, 
due to diffraction, after the shock passes the end of the obstacle. The orientation of the 
obstacle is not as critical as with the shear plates. The figure below illustrates these 
points 
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Figure 9 Flow field results for wedge obstacle 

The previous figure illustrates the effects of a wedge obstacle on the flow. The geometry 
of the wedge is such that one face is roughly perpendicular to the shock and a large 
pressure gradient builds on this face. This is shown by the closely spaced contour lines in 
the lower left comer between grid (4,4) and (7,2). The effects of this pressure gradient 
are that, the wedge absorbs a large amount of energy, the flow is split and a shadow or 
wake forms inside the faces of the wedge. The pressure gradient downstream of the 
wedge is ‘softened’, but its magnitude is, once again, similar to the unimpeded shock. 
Once again, following passage of the shock past the end of the wedge diffraction begins 
and the pressure begins to fill the shadow. This result is not unexpected. 

Since the orientation of the wedge is not as critical to its performance as that of the shear 
plates, wedges would be more a practical protection than the other obstacles investigated. 
However, due to the low viscosity of air the wedge would, in all likelihood have to be 
incorporated as an element on the exterior of the structure, and not a separate element. 
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This incoiporation would lead to a structure with wedge shaped faces which, it is already 
known to demonstrate good resistance to surface loading from many sources including 
blasts, waves, running water etc. Finally this form of protection is neither inexpensive 
nor aesthetically pleasing. 



CHAPTER 6: COST ANALYSIS 



The cost analyses for the obstacles investigated are moot. Since the pole obstacle does 
not work as hoped, unless many poles are placed closely together, in which case a wall 
would most likely be cheaper, I will present a synopsis of the available cost analyses for 
the conventional approach to blast protection. The following text is a compilation of 
others work and is not to be considered my own. These are my words but not my ideas or 
effort. 

Nuclear Disasters and the Built Environment contains an overview of the results of a 
testing completed in the 1950’s with above ground nuclear weapons. The results show 
that an overpressure of 5 psi (34.4 kPa) will completely demolish a conventional wood 
frame house and an oveipressure of 1.7 psi (1 1.7 kPa) will result in serious damage. A 
later test series included strengthened wood frame houses. The first house, which was 
subjected to 4 psi (27.6 kPa) oveipressure, sustained serious structural damage but 
remained standing. The second house was exposed to 2.6 psi (17.9 kPa) and suffered 
damage similar to the unreinforced house exposed to 1.7 psi. Also included in the second 
test were two brick houses. These were found to be of the same strength as the wood 
frame houses. One brick house was exposed to 5 psi (34.4 kPa) and sustained damage 
similar to the wood frame house. The second house received an oveipressure of 1 .7 psi 
(11.7 kPa) and sustained much less damage than the wood frame house. 

The strengthening of the houses for the test resulted in a 5% construction premium over 
unreinforced construction. 

The table below compiles the strengths of some common building materials. As can be 
seen the overpressure required for failure is quite low. 



44 



Table 1 1 Strength of common building materials 



Structural element 


Overpressure to cause failure 
(psi) 


Unreinforced glass 


0.5-1 .0 


Coirugated steel or aluminum siding 


1. 0-2.0 


Brick wall, unreinforced 


3.0-10.0 


Standard house construction 


1. 0-2.0 


Concrete wall panels 8”- 12” thick 


1. 5-5.5 



Protecting Buildings from Bomb Damage contains an extensive and thorough cost 
analysis of hardening a new building to resist exterior explosions. I will summarize the 
results of the analysis. The model was based upon a 5% construction premium for 
hardening the building and also assumed a minimum 10% return on investment. The 
construction premium is the additional cost that would be incurred if the blast hardening 
were to be included in the construction of a building. The construction premium 
assumption is in keeping with the result from the tests completed in the 1950’s. The 
output from the model was the lease premium based upon these assumptions. These 
results are shown in the table below. 



45 



Table 12 Constructi on-lease premium 



Assume construction premium 


Resulting lease premium 


3% 


2.69% 


5% 


3.46% 


7% 


4.23% 



The lease premium in the table above is the additional amount that would have to be 
charged lessees in a blast hardened building. This lease premium would recover the 
construction premium at the 10% return on investment. The table clearly shows that the 
resulting lease premium is less than the construction premium. Therefore it can be 
economically feasible to provide blast protection in new construction, if the developer or 
customer deems it necessary. 
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APPENDIX A: MACCORMACK 1-D CODE LISTING 



% 

% Some constants needed for calculations (UNITS) 

%Courant number 

CFL=input('Enter the Courant number(0<CFL<l 0.5 is a good starting point): '); 
if CFL>1. 

error('The Courant number must be less than 1 !') 
elseif CFL<=0. 

error('The Courant number must be larger than zero!') 
end 



%gamma for air 

g=l-4; 

gi=g-i; 

% 

%rho = initial density 

(kg/m A 3) 

rho=1.614; 

%Pin = initial pressure (Pa) 

Pin=101325; 

% 

%T=initial temperature 
(K) 

T 1=298; 

% 



%calculate the speed of sound at the initial pressure 
c=( 1 ,4*Pin/rho) A ( 1 / 2 ); 

% 



(m/s) 



%h= distance over which to perform calculations (m) 

h=input ('Enter the distance in meters over which to perform calculations: '); 
if h<=0. 



error('You must enter a nonzero distance!') 
end 



%stop= number of position steps per distance 

stop=input ('Enter the total number of x grid points(101 is a good starting point): '); 
if stop<=0. 
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enor('You must enter a nonzero number of x grid points!') 
end; 

% 

%calculate the x step size 
dx=h/(stop-l); 

% 

% Print output of values on screen 

% 

fprintf(' \n') 

fprintf(' Calculation values 

\n') 

fprintf('The Courant number is % 1 .2f\n',CFL) 

fprintf('The distance to calculate over is %3.2f m\n',h) 

fprintf('The number of x grid points is %4.0f\n',stop-l) 

fprintf('The x step size is %1.4f m\n',dx) 

fprintf(' \n') 

fprintf(’ ***************************************************************** 
*\n') 

fprintf(' Conditions ahead of the shock 

\n') 



fprintf('Gamma is %1 . 1 f\n',g) 

fprintf('The initial pressure is %3.3g Pa\n',Pin) 

fprintf('The initial density is %1.4fkg/m A 3\n',rho) 

fprintf('The initial temperature is %3.3g K\n',Tl) 

fprintf('The speed of sound is %3.4f m/s\n',c) 

fprintf( 1 **************************************************************** * 
*\n') 

% 

endt=input('Enter the number of time steps to calculate over: '); 
if endt<=0 



error('You must enter a number of time steps larger than zero!') 
end 

% 

%define vector size 



oldp=[l:stop]*0.0; 
newp=[ 1 :stop]*0.0; 
oldr=[l :stop]*0.0; 
newr=[l :stop]*0.0; 
barr=[l:stop]*0.0; 
oldru=[l:stop]*0.0; 
newru=[l:stop]*0.0; 
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bami=[l: stop] *0.0; 
newu=[l: stop] *0.0; 
oldE=[l :stop]*0.0; 
newE=[ 1 :stop]*0.0; 
barE=[l :stop]*0.0; 
loc=[l:stop]*0.0; 

% 

% Convert x grid points to actual distances 



for x=l :stop 

loc(x)=(x-l)*dx; 

end 

% 

% Calculate Peak staic overpressure Ps based on charge size 

% and range to charge the driving force for the 

% calculations to follow 

% 

% i - range to charge 

(m) 

r=input('Enter the range to the charge in meters: '); 
if r<=0 

error('You must enter a range larger than zero!') 
end 

% Print out table of explosive types 



% T=type of explosive 

fprintf(' \n') 

fprintf(' Choose the explosive using the table below\t\n') 

fprintf('Compound B l\n') 

fprintf('RDX (Cyclonite) 2\n') 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin (liquid) 4\n') 

fprintf('TNT 5\n') 

fprintf('Blasting Gelatin 6\n') 

fprintf(’60 percent Nitroglycerin dynamite 7\n') 

fprintf('Semtex 8\n') 

fprintfC — — \n') 

% 



% Ask for whiat type of explosive to use 

T=input('enter the type of explosive to use: '); 
if T=1 



S=1 . 148; 
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elseif T==2 



el seif T==3 



elseif T=4 



elseif T==5 



elseif T==6 



elseif T~7 



elseif T==8 



S= 1.185; 
S=1 .256; 
S=1 .481; 
S=1.000; 
S=1.000; 
S=0.600; 



S=1 .250; 

else 

eiTor('You must enter an explosive type!') 
end 

fprintf('The TNT conversion factor for this explosive is: %1 ,3f\n',S) 
% 



% Ask what charge mass to use (kg of TNT) 

% 

M=input(' Enter the mass of the charge in kilograms: '); 
if M<0 

error('You must enter a mass larger than zero!') 
end 



W=M*S; 

% 

% calculate peak static overpressure (Bar) 



z=r/(W A (l/3)); 

Ps=(6.7/(z A 3)), 

if Ps< 1 0 

Ps=((0.975/z)+(1.455/(z A 2))+(5.85/(z A 3)))-0.019; 
else Ps=(6.7/(z A 3)); 
end 



if Ps<0. 1 
end 



Ps=Ps; 



%- 

% 

% 



Convert Ps(the static overpressure) from Bar (Pa) 
to Pascals from overpressure to absolute pressure 
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% and apply conversion factor for ground burst 

% 

Ps= 1 .8 *((Ps* 1 00000)+Pin); 



% 

% Print output of values on screen 



fprintf(' \n’) 

fprintf(' Explosive values 

\n') 

fprintf('The range to the charge from the x=0 point is %4.2f m\n',r) 

fprintf('The mass of the charge in TNT equivalents is %4.3f kg\n',W) 

fprintf('The static overpressure of the charge is %1.3g Pa\n',Ps-Pin) 

fprintf(' \n') 



% Calculate intial velocity using Ps 

% 

Mach=sqrt(((Ps/Pin)- 1 )*((g+ 1 )/(2 *g))+ 1 ); 

Vin=c*(2/(g+ 1 )*(Mach-( 1/Mach))); 

% Print output of values on screen 

fprintf('The initial wave velocity is %4.4f m/s\n',Vin) 

fprintf('The shock wave Mach number is %2.4f\n',Mach) 

fprintf(' \n') 

% 

% Calculate values behind the shock 

% 

rho2=rho/( 1 -(2/(g+ 1))*(1-(1 /Mach A 2))); 

T2=Ps/(rho2*(83 14.5 1/29)); 
c2=(c A 2 *T2/T 1 ) A 0. 5; 

fpri ntf(’ ***************************************************************** 
*\n') 

fprintf(' Conditions behind the shock 

\ri) 

fprintf('The pressure is %1.3g Pa\n',Ps) 

fprintf('The density is %1.4fkg/m A 3\n',rho2) 

fprintf('The temperature is %3.3g K\n',T2) 

fprintf('The speed of sound is %3.4f m/s\n',c2) 

fprin tf( ****************************************************************** 
*\n’) 



% 

% 

% 



Define initial conditions 
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for x= Estop 

oldp(x)=Pin; 

oldr(x)=rho; 

oldru(x)=0; 

oldE(x)=Pin/gl; 

end 

% 

% Define boundary conditions at grid point x= 1 



oldp(l)=Ps; 

oldr( 1 )=rho/( 1 -(2/(g+ 1 ))*( 1 -( 1 /Mach A 2))); 
oldru( 1 )=oldr( 1 )*' Vin; 
oldE(l)=(Ps/gl)+(0.5*oldru(l) A 2/oldr(l)); 



% 

% Calculate the new speed of sound based on new pressure 

% and density values 

% 

newc=sqrt( 1 ,4*oldp(x)/oldr(x)); 

% 

% Calculate the initial time step 

% 

dt=CFL*dx/(newc+oldru( 1 )/oldr( 1 )); 

fprintf(’The inital time step is %1.5e s\n',dt) 

% 

% Initialize total time counter 

% 

tt=0; 

% 

% Define predictor values at x equal 1 

% 

barr( 1 )=oldr( 1 )-(dt)*(oldm(2)-oldm( 1 )); 
barru(l)=oldru(l)-(dt)*(((oldm(2) A 2/oldr(2))-... 

(oldru( 1 ) A 2/oldr( 1 )))+(oldp(2)-oldp( 1 ))); 
barE(l)=oldE(l)-(dt)*((oldE(2)+oldp(2))*oldru(2)/... 
oldr(2)-(oldE( 1 )+oldp( 1 ))*oldru( 1 )/oldr( 1 ) ); 



% 

% Adjust x grid end point for finite differencing 



stop=stop-l; 
% 



54 



% Finite difference for density, momentum and energy 

% 



for t=l:endt 

dtmin=l.; 
CFLmin=CFL; 
dtdx=dt/(dx); 
%Calculate total time 
tt=tt+dt; 
for x=2:stop 



% 

% 

% 

% 

% 

%• 

% 



% — 
% 



% 

% 



% — 
% 

% — 
% 



% 



PREDICTOR 

Solve as a function of time and position for increasing 
time 

Density(r) 

ban-(x)=oldr(x)-(dtdx)*(oldru(x+ 1 )-oldru(x)); 
Momentum(ru) 

bami(x)=oldru(x)-(dtdx)*((oldru(x+l) A 2/oldr(x+l ))-... 
(oldru(x) A 2/oldr(x))+(oldp(x+l)-oldp(x))); 

Energy(E) 

barE(x)=oldE(x)-(dtdx)*((oldE(x+ 1 )+oldp(x+ 1 ))*. .. 
(oldru(x+ 1 )/ol dr(x+ 1 ))-(oldE(x)+ol dp(x)) . . 

(ol dru( x)/ol dr(x))); 



CORRECTOR 



Density(r) 

newr(x)=(barr(x)+oldr(x))/2-(dtdx/2)*(bairu(x)-ban'u(x-l )); 

Left wall boundaiy conditions 

newr(stop+ 1 )=newr(stop- 1 ); 
newr(stop)=newr(stop- 1 ); 
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% 

% Momentum(ru) 

newru(x)=(barru(x)+oldru(x))/2-(dtdx/2)*((banu(x) A 2/... 
ban(x))-(barru(x- 1 ) A 2/barr(x-l))+(g 1 *(barE(x)-... 

(. 5 *(barru(x) A 2/barr(x)))))-(g 1 *(barE(x- 1 )-( . 5 * . . . 
(barru(x- 1 ) A 2/barr(x- 1 )))))); 

% Left wall boundary conditions 

newru(stop)=0; 

newru(stop+l)=0; 



% 

% Energy(e) 

newE(x)=(barE(x)+oldE(x))/2-(dtdx/2)*((barru(x)/... 
barr(x))*(barE(x)+(g 1 *(barE(x)-(.5 *(barru(x) A 2/. .. 
baiT(x))))))-((bairu(x- 1 )/barr(x- 1 ))*(barE(x- 1 )+. . . 
(g 1 *(barE(x- 1 )-(. 5 *(banoi(x-l ) A 2/barr(x- 1 )))))))); 



% 

% Calculate pressure using the continuity equation for the 

% ideal gas law and the new values for energy, momentum 

% and density 

newp(x)=gl*(newE(x)-(.5*(newru(x) A 2/newr(x)))); 

newp(stop+ 1 )=newp(stop- 1 ) ; 



% 

% Calculate velocity using the new values for momentum 

% and density 

newu(x)=newru(x)/newr(x); 

% 

% Check stability 

% Calculate the new speed of sound 



newc=sqrt(g*newp(x)/newr(x)); 
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% 

% Calculate the new time step dt 



dtnew=dx/ (newu(x)+newc); 
if dtnew < dtmin 

dtmin=dtnew; 

end 

CFLnew=(dtmin/dx)*((newc+oldru(x)/oldr(x))); 
if CFLnew<CFLmin 

CFLmin=CFLnew; 

end 

% 

% Return for next x 

end 



% Apply numerical boundaiy scheme to calculate all values 

% atx=l andx=stop+l 

newr(l)=oldr(l); 
newE(l)=oldE( 1); 
newru( 1 )=oldru( 1 ); 
newu( 1 )=oldru( 1 )/oldr( 1); 
newp(l)=oldp(l); 

barr( 1 )=oldr( 1 )-(dt)*(oldru(2)-oldru( 1 )); 
barru( 1 )=oldru( 1 )-(dt)*((oldru(2) A 2/ oldr(2)-. . . 

oldru( l) A 2/oldr( l))+(oldp(2)-oldp(2))); 
barE( 1 )=oldE( 1 )-(dt)*((oldE(2)+oldp(2))*oldru(2)/... 

oldr(2) -(oldE(l)+oldp( l))*oldm( l)/oldr(l) ); 
newE(stop+ l)=oldE(stop+ 1 ); 
newu(stop+ 1 )=newru(stop+ 1 )/newr(stop+ 1 ); 



% Update time step and CFL(Courant number) 



dt=dtmin; 

CFL=CFLmin; 

fprintf('CFL= %1 .3g\tdelta t= %1.5e s\tTotal time= %1.5e s\n',CFL,dt,tt) 



% — 
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% Output the values for pressure, density and velocity 

% using graphs 

% 

% Pressure 

subplot(3,l,l), plot(loc,newp), ylabel('Absolute Pressure (Pa)') 

% 

% Density 

subplot(3,l,2), plot(loc,newr), ylabel('Density (kg/m A 3)') 

% Velocity 

subplot(3,l,3), plot(loc,newu), ylabel('Velocity (m/s)'),... 
xlabel('Distance (m)') 



% 

pause(.OOOl) 

% 

%Replace old values with new values and begin next time step 

oldp=newp; 

oldt-newr; 

oldru=newru; 

oldE=newE; 

% 

% Return for next t 

end 



APPENDIX B: MACCORMACK QUASI 2-D CODE LISTING 

clear all 



% 



Some constants needed for calculations (UNITS) 



%Courant number 

CFL=input('Enter the Courant number(0<CFL<l 0.5 is a good starting point): ') 
if CFL>1. 

error('The Courant number must be less than 1 !') 
elseif CFL<=0. 

error('The courant number must be larger than zero!') 
end 

%gamma for air 

g=1.4; 

gl=g-l; 



% 

%rho = initial density 

(kg/m A 3) 

rho=1.614; 



% 

%Pin = initial pressure (Pa) 

Pin=101325; 



% 

%T=initial temperature 

(K) 

T 1=298; 



% 

%calculate the speed of sound at the initial pressure (m/s) 

c=(1.4*Pin/rho) A (l/2); 



% 

%h= distance over which to perform calculations (m) 

h=input ('Enter the distance in meters over which to perform calculations: '); 
if h<=0. 

error('You must enter a nonzero distance!') 
end 
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% 

%stop= number of position steps per distance 

stop=input ('Enter the total number of grid points in the x,y directions(3 1 is a good 
starting point): '); 

if stop<=0. 

error('You must enter a nonzero number of x gnd points!') 
end; 



% 

%calculate the x and y step size 

dx=h/(stop-l); 

dy=h/(stop-l); 

% 



% Print output of values on screen 

fprintfC \n') 

fprintf(' Calculation values 

\n') 

fprintf('The Courant number is %1 ,2f\n',CFL) 

fprintf('The x distance to calculate over is %3.2f m\n',h) 

fprintf('The y distance to calculate over is %3.2f m\n',h) 

fprintf('The number of x grid points is %4.0ftn',stop-l) 

fprintf('The number of y grid points is %4.0ftn',stop-l) 

fprintf('The x step size is %1 .4f m\n’,dx) 

fprintf('The y step size is %1 .4f m\n',dx) 

fprintf(' \n') 

*\n') 

fprintf(' Conditions ahead of the shock 

\n') 

fprintf('Gamma is %1. lf\n',g) 

fprintf('The initial pressure is %3.3g Pa\n',Pin) 

fprintf('The initial density is %1.4f kg/m A 3\n',rho) 

fprintf('The initial temperature is %3.3g K\n',Tl) 

fprintf('The speed of sound is %3.4f m/s\n',c) 

*\n) 



endt=input('Enter the number of time steps to calculate over: '); 
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if endt<=0 

error('You must enter a number of time steps larger than zero!') 
end 



% 

% define vector size 



for x=l:stop 

for y= Estop 

op(x,y)=0.0; 

np(x,y)=0.0; 

Olr(x,y)=0.0; 

nr(x,y)=0.0; 

br(x,y)=0.0; 

Olru(x,y)=0.0; 

nru(x,y)=0.0; 

Olrv(x,y)=0.0; 

nrv(x,y)=0.0; 

bru(x,y)=0.0; 

nu(x,y)=0.0; 

brv(x,y)=0.0; 

nv(x,y)=0.0; 

oE(x,y)=0.0; 

nE(x,y)=0.0; 

bE(x,y)=0.0; 



end 

end 

% 

% 

% 

% 

o / 0 

% 



% 



Calculate Peak static overpressure Ps based on charge 
size and range to charge the driving force for the 
calculations to follow 



r=range to charge 

(m) 



r=input('Enter the range to the charge in meters: '); 
if r<=0 

error(’You must enter a range larger than zero!') 
end 
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% Print out table of explosive types 

% T=type of explosive 



fprintf(' \n') 

fprintf(' Choose the explosive using the table below\t\n') 

fprintf('Compound B l\n') 

fprintf('RDX (Cyclonite) 2\n') 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin (liquid) 4\n') 

fprintf('TNT 5\n') 

fprintf('Blasting Gelatin 6\n') 

fprintf('60 percent Nitroglycerin dynamite 7\n') 

fprintf('Semtex 8\n') 

fprintf(' -\n') 



% Ask for what type of explosive to use 

% 



T=input('enter the type of explosive to use: '); 
if T=1 



el seif T=2 
el seif T=3 
elseif T=4 
elseif T=5 
elseif T==6 
elseif T==7 
elseif T=8 



S=l. 148; 
S=l. 185; 
S=1.256; 
S= 1 .48 1 ; 
S= 1.000; 
S=1.000; 
S=0.600; 
S=1.250; 



else 

error('You must enter an explosive type!') 
end 

fprintf('The TNT conversion factor for this explosive is: %1.3ftn',S) 



% 



% 



62 



Ask for mass of charge 
(kg) 



M=input(' Enter the mass of the charge in kilograms: '); 
if M<0 

error('You must enter a mass larger than zero!') 
end 

W=M*S; 

% 

% Calculate peak static overpressure (B) 

% 



z=r/(W A (l/3)); 

Ps=(6.7/(z A 3)); 

if Ps< 1 0 

Ps=((0.975/z)+(1.455/(z A 2))+(5.85/(z A 3)))-0.019; 
else Ps=(6.7/(z A 3)); 
end 



if Ps<0. 1 
end 



Ps=Ps; 



% 

% Convert Ps(the static oveipressure) from Bar (Pa) 

% to Pascals from overpressure to absolute 

% pressure and apply conversion factor for 

% ground burst 

% 

Ps=l ,8*((Ps* 1 00000)+Pin); 



% 

% Print output of values on screen 

% 



fprintf(' Explosive values 

\n') 

fprintf('The range to the charge from the x=0 point is %4.2f m\n',r) 
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fprintf('The mass of the charge in TNT equivalents is %4.3f kg\n',W) 

fprintf('The static oveipressure of the charge is %1.3g Pa\n',Ps-Pin) 

fprintf(' \n') 



% Calculate intial velocity using Ps 

% 



Mach=sqrt(((Ps/Pin)-l)*((g+l)/(2*g))+l); 
Vin=c*(2/(g+ 1 )*(Mach-( 1 /Mach))); 

% 

fprintf('The initial wave velocity is 

fprintf('The shock wave Mach number is... 
fprintf(' 



Print output of values on screen 

%4.4f m/s\n',Vin) 

%2.4f\n',Mach) 

.\n') 



% Calculate values behind the shock 



rho2=rho/( 1 -(2/(g+ 1 ))*( 1 -( 1 /Mach A 2))); 

T2=Ps/(rho2*(83 1 4.5 1/29)); 
c2=(c A 2 *T2/T 1 ) A 0 . 5; 

fprintfO* **************************************************************** 
*\n') 

fprintf(' Conditions behind the shock 

\n') 

fprintf('The pressure is %1.3g Pa\n',Ps) 

fprintf('The density is %1 ,4f kg/m A 3\n',rho2) 

fprintf('The temperature is %3.3g K\n',T2) 

fprintf('The speed of sound is %3.4f m/s\n',c2) 

fpri ntf(* ***************************************************************** 
*\n') 

% 



% 




% 




% 


— 


for 


x=l :stop 


for 


y=l:stop 

op(x,y)=Pin; 

01r(x,y)=rho; 

Olru(x,y)=0; 



Define initial conditions 
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Olrv(x,y)=0; 

oE(x,y)=Pin/gl; 

end 

end 

% Define boundary conditions at grid point x=l y=l 

% 

y=i; 

for x= Estop 

op(x,y)=Ps; 

01r(x,y)=rho/( 1 -(2/(g+ 1 ))*( 1 -( 1 /M ach A 2))); 
Olru(x,y)=0; 

01rv(x,y)=01r(x,y)*Vin; 

oE(x,y)=(Ps/gl)+(0.5*(01ru(x,y) A 2+01rv(x,y) A 2)/... 

01r(x,y)); 



end 



% 

% Calculate the new speed of sound based on new pressure 

% and density values 

nc=sqrt(g*op(x,y)/01r(x,y)); 

% Calculate the initial time step 

% 



dt=CFL*( 1 /(abs(01ru( 1 , l)/01r(l , 1 ))/dx+abs(01rv( l,l)/01r( 1 , 1 ))/... 
dy+nc*(l/dx A 2+l/dy A 2) A (0.5))); 



fprintf('The inital time step is %1.5e s\n',dt) 

% 



tt=0; 



Initialize total time counter 
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% 

% Finite difference for density, momentum and energy 

% 

% Adjust x and y grid end point for finite differencing 

% 



stop=stop- 1 ; 



for t= 1 : endt 

% 

% NBS for column 1 and row 1 



dtmin=l.; 

dtdx=dt/dx; 

dtdy=dt/dy, 



% 

% Calculate total time 



tt=tt+dt; 



for x=l:stop 
for y=l:stop 






% PREDICTOR 

% — 

% 

% 



Solve as a function of time and position for 
increasing time 
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Density(r) 



% 

% 



br(x,y)=01r(x,y)-... 

(dt*(((01ru(x+ 1 ,y)-01ru(x,y))/dx)+. . . 

((01rv(x,y+l)-01rv(x,y))/dy))); 

br(x,stop+ 1 )=br(x,stop); 
br(stop+ 1 ,y)=br(stop,y); 
br(stop+ 1 ,stop+ 1 )=br(stop,stop); 

% 

% 

% 



bru(x,y)=01ru(x,y)-dt*. . . 

((1/dx*... 

(01ru(x+l ,y) A 2/01r(x+l ,y)-... 
01ru(x,y) A 2 /01r(x,y) +... 

(op(x+l ,y)-op(x,y) )))+... 

(1/dy*... 

(01rv(x,y+ 1 )*01ru(x,y+ 1 )/01r(x,y+ 1 )-. . . 
01rv(x,y) *01ru(x,y) /01r(x,y) ))); 

bru(x,stop+ 1 )=bm(x,stop); 
bru(stop+l ,y)=bm(stop,y); 
bru(stop+ 1 ,stop+ 1 )=bru(stop,stop); 



% 

% 



brv(x,y)=01rv(x,y)-dt*... 

((1/dx*... 

(01ru(x+l ,y)*01rv(x+l ,y)/01r(x+l ,y)-... 
01ru(x,y) *01rv(x,y) /01r(x,y) ))+... 
(1/dy*... 

(01rv(x,y+ 1 ) A 2/0h'(x,y+ 1 )-. . . 
01rv(x,y) A 2 /01r(x,y)+... 
(op(x,y+l)-op(x,y) )))); 



Momentum(ru) 



Momentum(rv) 



brv(stop+l ,y)=brv(stop,y); 
brv(x,stop+ 1 )=brv(x,stop); 
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brv(stop+ 1 ,stop+ 1 )=brv(stop,stop); 

% 

% Energy(e) 



bE(x,y)=oE(x,y)-dt*... 

((1/dx*... 

((oE(x+ 1 ,y)+op(x+ 1 ,y)) *01ru(x+ 1 ,y)/01r(x+ 1 ,y) 
(°E(x,y)+op(x,y) )*01ru(x,y) /01r(x,y) )+... 
1/dy*... 

((oE(x,y+ 1 )+op(x,y+ 1 )) *01rv(x,y+ 1 )/01r(x,y+ 1 ) 
(oE(x,y)+op(x,y) )*01rv(x ,y) /01r(x,y) ))); 

bE(stop+ 1 ,y)=bE(stop,y); 
bE(x,stop+ 1 )=bE(x,stop); 
bE(stop+ 1 ,stop+ 1 )=bE(stop,stop); 



% return for next x 
end 

% return for next y 
end 



% 

nr(:,l)=01r(:,l); 

nru(:,l)=01ru(:,l); 

nrv(:,l)=01rv(:,l); 

nE(:,l)=oE(:,l); 

for x=2:stop+l 

for y=2:stop+l 

% CORRECTOR 

% Density(r) 



nr (x,y)=(br(x,y)+01r(x,y))/2-... 

( dtdx/2*(bru(x,y)-bru(x-l,y))+... 
dtdy/2*(brv(x,y)-brv(x,y- 1 ))); 
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nr( 1 ,y)=nr(2,y); 

Left wall boundary conditions 



% 



% 



% Momentum(ru) 

% 

Px=gl*(bE(x,y)-0.5*((bru(x,y) A 2/br(x,y))+(brv(x,y) A 2/br(x,y)))); 

Pxm=g 1 *(bE(x- 1 ,y)-0. 5 *((bru(x- 1 ,y) A 2/br(x- 1 ,y))+(brv(x- 1 ,y) A 2/br(x- 1 ,y)))); 

nru(x,y)=(bru(x,y)+01ru(x,y))/2-... 

((dtdx/2*... 

((bru(x,y) *bru(x,y) /br(x,y) +Px)-... 

( bru(x- 1 ,y)*bru(x-l ,y)/br(x- 1 ,y)+Pxm)))+... 

(dtdy/2*... 

( bru(x,y) *brv(x,y) /br(x,y)-... 
bru(x,y- 1 )*bi"v(x,y- 1 )/br(x,y- 1 ) ))) ; 



% 



nru( 1 ,y)=nru(2,y); 



% Momentum(rv) 

0/ 0 



Py=gl*(bE(x,y)-0.5*((bru(x,y) A 2/br(x,y))+(bi-v(x,y) A 2/br(x,y)))); 

Pym=g 1 *(bE(x,y- 1 )-0. 5*((bm(x,y- 1 ) A 2/br(x,y- 1 ))+(brv(x,y- 1 ) A 2/br(x,y- 1 )))); 

nrv(x,y)=(brv(x,y)+01rv(x,y))/2-... 

((dtdx/2*... 

( bru(x,y) *bi"v(x,y) /br(x,y) 
bm(x-l,y)*brv(x-l ,y)/br(x-l,y) ))+... 

(dtdy/2*... 

( (brv(x,y)*brv(x,y)/br(x,y)+ Py )-... 

(brv(x,y- 1 ) *biw(x,y- 1 )/br(x,y- 1 )+Pym)))), 






nrv( 1 ,y)=nrv(2,y); 
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% Left wall boundary conditions 



% 



% Energy(E) 



nE(x,y)=(bE(x,y)+oE(x,y))/2-... 

((dtdx/2*... 

((bru(x,y) /br(x,y) *(bE(x,y) +Px )) 

(bru(x- 1 ,y)/br(x-l ,y)*(bE(x-l ,y)+Pxm)))) +... 
(dtdy/2*... 

((brv(x,y) /br(x,y) *(bE(x,y)+Py )) 

(brv(x,y- 1 )/br(x,y- 1 ) *(bE(x,y- 1 )+Pym)))) ); 

% 

% 

% Right wall boundary conditions 



nr(x,stop+ 1 )=nr(x,stop- 1 ); 

% nr(x,stop)=nr(x,stop- 1 ); 

nrv(x,stop)=0; 
nrv(x,stop+l)=0; 
nE(x,stop+l)=oE(x,stop+l); 

nE(l,y)=nE(2,y); 

% 

%retum for next x 
end 

%retum for next y 
end 

% 



for 
for 
% — 
% 

% 

% 

% — 



y=l:stop+l 

x=l:stop+l 



Calculate pressure using the continuity equation for 
the ideal gas law and the new values for energy, 
momentum and density 
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np(x,y)=(g 1 ) *(nE(x,y)-(. 5 *((nru(x,y) A 2/. . . 
m(x,y))+(m-v(x,y) A 2/nr(x,y))))); 

np(x,stop+l)=np(x,stop-l); 



% 

% Calculate velocity using the new values for momentum 

% and density 

% 



nu(x ,y)=nru(x,y)/ nr(x,y) ; 
nv(x,y)=nrv(x,y)/nr(x,y); 



% 

% Check stability 

% Calculate the new speed of sound 

% 



nc=sqrt(g*np(x,y)/nr(x,y)); 

% 

% 

% 

0/ 0 



Calculate the new time step dt 



dtn=CFL*(l/(abs(nu(x,y))/dx+abs(nv(x,y))/dy+... 

nc*(l /dx A 2+ 1 /dy A 2) A (0. 5))); 
if dtn < dtmin 

dtmin=dtn; 

% end dtmin if statement 
end 

0/ 0 

% 



% 

%retum for next x 
end 

o / 0 

%retum for next y 
end 
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% 

% Update time step and CFL(Courant number) 

% 



dt=dtmin; 

fprintf('delta t= %1.5e s\tTotal time= %1.5e s\n',dt,tt) 



% 

% Output the values for pressure, density and velocity 

% using graphs 

% Pressure 



colormap(cool) 
surf(np) 
hold on 

pcolor((np)-Pin) 
title('Absolute Pressure') 
xlabel('X grid points') 
ylabel('Y grid points') 
zlabel('Pascals') 

pause(.Ol) 
hold off 

% 

% Replace old values with new values and begin 

% next time step 

% 



op=np; 

01r=nr, 

01ru=nm; 

01rv=nrv; 

oE=nE; 

% 



% Return for next t 

end 



APPENDIX C: MACCORMACK 2-D CODE LISTING 

clear all 

% 



% 



Some constants needed for calculations (UNITS) 



%Courant number 

CFL=input('Enter the Courant number(0<CFL<l 0.5 is a good starting point): '); 
if CFL>1. 

error(’The Courant number must be less than 1 !') 
elseif CFL<=0. 

error(’The courant number must be larger than zero!') 
end 

%gamma for air 

g=1.4; 

gl=g-1; 



% 

%rho = initial density 

(kg/m A 3) 

rho=1.614; 



% 

%Pin = initial pressure (Pa) 

Pin=l 0 1325; 

%T=initial temperature 

(K) 

T 1=298; 



% 

%calculate the speed of sound at the initial pressure (m/s) 

c=( 1 ,4*Pin/rho) A (l/2); 



% 

%h= distance over which to perform calculations (m) 

h=input (’Enter the distance in meters over which to perform calculations: ’); 
if h<=0. 

enor('You must enter a nonzero distance!') 
end 
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% 

%stop= number of position steps per distance 

stop=input ('Enter the total number of grid points in the x,y directions(3 1 is a good 
starting point): '); 

if stop<=0. 

error('You must enter a nonzero number of x grid points!') 
end; 



%calculate the x and y step size 

dx=h/(stop-l); 

dy=h/(stop-l); 



% Print output of values on screen 

fprintf(' -W) 

fprintf(' Calculation values 

\n') 

fprintf('The Courant number is %1.2f\n',CFL) 

fprintf(’The x distance to calculate over is %3.2f m\n',h) 

fprintf('The y distance to calculate over is %3.2f m\n',h) 

fprintf('The number of x grid points is %4.0f\n', stop-1) 

fprintf('The number of y grid points is %4.0f\n',stop-l) 

fprintf('The x step size is % 1 4f m\n',dx) 

fprintf('The y step size is %1.4f m\n',dx) 

fprintf(' \n') 

fprintf(’ ***************************************************************** 
*\n') 

fprintf(' Conditions ahead of the shock 

\n') 



fprintf('Gamma is %l.lf\n',g) 

fprintf('The initial pressure is %3.3g Pa\n',Pin) 

fprintf('The initial density is %1.4fkg/m A 3\n',rho) 

fprintf('The initial temperature is %3.3g K\n',Tl) 

fprintf('The speed of sound is %3.4f m/s\n',c) 

fprintf(' ***************************************************************** 
*\n') 



% — 
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endt=input('Enter the number of time steps to calculate over: '); 
if endt<=0 

eiTor('You must enter a number of time steps larger than zero!') 
end 



% — 
% 



define vector size 



for x= Estop 

for y=l:stop 

op(x,y)=0.0; 

np(x,y)=0.0; 

Olr(x,y)=0.0; 

nr(x,y)=0.0; 

br(x,y)=0.0; 

Olru(x,y)=0.0; 

nru(x,y)=0.0; 

Olrv(x,y)=0.0; 

nrv(x,y)=0.0; 

bru(x,y)=0.0; 

nu(x,y)=0.0; 

brv(x,y)=0.0; 

nv (x,y)=0.0; 

oE(x,y)=0.0; 

nE(x,y)=0.0; 

bE(x,y)=0.0; 



end 

end 



% 


— 


% 


Calculate Peak static overpressure Ps based on charge 


% 


size and range to charge the driving force for the 


% 


calculations to follow 


% 


— 


% 


r=range to charge 




(m) 


% 





1 -input('Enter the range to the charge in meters: '); 
if r<=0 

eiTor('You must enter a range larger than zero!') 



76 



end 

% Print out table of explosive types 

% T=type of explosive 

fprintf(' \n') 

fprintf(' Choose the explosive using the table below\t\n') 

fprintf('Compound B l\n') 

fprintf('RDX (Cyclonite) 2\n') 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin (liquid) 4\n') 

fprintf('TNT 5\n') 

fprintf('Blasting Gelatin 6\n') 

fprintf('60 percent Nitroglycerin dynamite 7\n') 

fprintf('Semtex 8\n') 

fprintf(' \n') 

% 

% Ask for what type of explosive to use 

% 

T=input('enter the type of explosive to use: '); 
if T=1 

S=1.148; 

elseif T==2 

S= 1.185; 

elseif T==3 

S=1.256; 

elseif T==4 

S= 1 .48 1 ; 

elseif T==5 

S= 1.000; 

elseif T=6 

S= 1.000; 

elseif T==7 

S=0.600; 

elseif T==8 

S=1.250; 

else 

error('You must enter an explosive type!') 
end 

fprintf('The TNT conversion factor for this explosive is: %1.3f\n',S) 



% — 



% 
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Ask for mass of charge 
(kg) 



M=input('Enter the mass of the charge in kilograms: '); 
if M<0 

error('You must enter a mass larger than zero!') 
end 

W=M*S; 

% Calculate peak static overpressure (B) 

% 



z=r/(W A (l/3)); 

Ps=(6.7/(z A 3)); 

if Ps<10 

Ps=((0.975/z)+(1.455/(z A 2))+(5.85/(z A 3)))-0.019; 
else Ps=(6.7/(z A 3)); 
end 



if Ps<0. 1 
end 



Ps=Ps; 



% 

% 

% 

% 

% 

% 



Convert Ps(the static overpressure) from Bar (Pa) 
to Pascals from overpressure to absolute 
pressureand apply conversion factor for 
ground burst 



Ps= 1 .8*((Ps*100000)+Pin); 



% 

% Print output of values on screen 



fprintf(' \n') 

fprintf(' Explosive values 

\n') 

fprintf('The range to the charge from the x=0 point is %4.2f m\n',r) 
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fprintf('The mass of the charge in TNT equivalents is %4.3f kg\n',W) 

fprintf('The static overpressure of the charge is %1.3g Pa\n',Ps-Pin) 



% Calculate intial velocity using Ps 

% 



Mach=sqrt(((Ps/Pin)- 1 )*((g+ 1 )/(2*g))+ 1 ); 
Vin=c*(2/(g+l)*(Mach-( 1/Mach))); 

% 

fprintf('The initial wave velocity is 

fprintf('The shock wave Mach number is... 
fprintf(' 



Print output of values on screen 

%4.4f m/s\n',Vin) 

%2.4f\n',Mach) 

\n') 



% Calculate values behind the shock 



rho2=rho/( 1 -(2/(g+ 1))*(1 -( l/Mach A 2))); 

T2=Ps/(rho2*(83 14.5 1/29)); 
c2=(c A 2 *T2/T 1 ) A 0. 5; 

fprintf(* **************************************************** ************* 
*\n') 

fprintf(' Conditions behind the shock 

W) 

fprintf('The pressure is %1.3g Pa\n',Ps) 

fprintf('The density is %1.4fkg/m A 3\n',rho2) 

fprintf('The temperature is %3.3g K\n',T2) 

fprintf('The speed of sound is %3.4f m/s\n',c2) 

fprintf(' ***************************************************************** 
*\n') 



% 

% Define initial conditions 

for x=l:stop 

for y=l:stop 

op(x,y)=Pin; 

01r(x,y)=rho; 

Olru(x,y)=0; 



79 



Olrv(x,y)=0; 

°E(x,y)=Pin/gl; 

end 

end 

% Define boundary conditions at grid point x=l y=l 

y=2; 

x=2; 

op(x,y)=Ps; 

01r(x,y)=rho/( l-(2/(g+ 1))*( 1 -( l/Mach A 2))); 
01ru(x,y)=01r(x,y)*Vin; 

01rv(x,y)=01r(x,y)*Vin; 

o E(x ,y)=(Ps/ g 1 )+(0. 5 *(Olru(x ,y) A 2+01rv(x,y) A 2)/. . . 
01r(x,y)); 



% 

% 

% 

% 



Calculate the new speed of sound based on new pressure 
and density values 



n c=s qrt(g *op(x ,y)/01r (x ,y)) ; 

% 

% Calculate the initial time step 

% 

dt=CFL*( 1 /(abs(01ru( 1 , 1 )/01r( 1 , 1 ))/dx+abs(01rv( 1 , 1 )/01r( 1 , 1 ))/... 
dy+nc*(l/dx A 2+l/dy A 2) A (0.5))); 



fprintf('The inital time step is 

% 

% 



,...%1 ,5e s\n',dt) 

Initialize total time counter 



tt=0; 
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% 

% Finite difference for density, momentum and energy 

% 

% Adjust x and y grid end point for finite differencing 

% 



stop=stop-l; 



for t=l:endt 

% NBS for column 1 and row 1 

% 



dtmin=l.; 

dtdx=dt/dx; 

dtdy=dt/dy; 



% 

% Calculate total time 

% 



tt=tt+dt; 



% 

for x=l :stop 
for y=l:stop 






% PREDICTOR 

% 

% Solve as a function of time and position for 

% increasing time 
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% 

% 



br(x,y)=01r(x,y)-... 

(dt*(((01m(x+l,y)-01ru(x,y))/dx)+... 

((01rv(x,y+l)-01rv(x,y))/dy))); 

br(x,stop+ 1 )=br(x,stop); 
br(stop+ 1 ,y)=br(stop,y); 
br(stop+ 1 ,stop+ 1 )=br(stop,stop); 

% 

% 

% 



bru(x,y)=01ru(x,y)-dt*. . . 

((1/dx*... 

(01ru(x+ 1 ,y) A 2/01r(x+ 1 ,y)-. . . 
01ru(x,y) A 2 /01r(x,y) +... 

(op(x+ 1 ,y)-op(x,y) )))+... 

(1/dy*... 

(01rv(x ,y+ 1 ) * 01ru(x,y+ 1 )/01r(x,y+ 1 )-. . . 
01rv(x,y) *01ru(x,y) /01r(x,y) ))); 

bru(x,stop+l)=bru(x,stop); 
bru(stop+ 1 ,y)=bm(stop,y); 
bru(stop+ 1 ,stop+l )=bru(stop,stop); 



% 

% 

% 



brv(x,y)=01rv(x,y)-dt* . . . 

((l/dx*... 

(01ru(x+ 1 ,y)*01rv(x+ 1 ,y)/01r(x+ 1 ,y)-. . . 
01ru(x,y) *01rv(x,y) /01r(x,y) ))+... 
(1/dy*... 

(01rv(x,y+ 1 ) A 2/01r(x,y+ 1 )-. . . 
01rv(x,y) A 2 /01r(x,y)+... 
(op(x,y+l)-op(x,y) )))); 



Density(r) 



Momentum(ru) 



Momentum(rv) 



brv(stop+ 1 ,y)=brv(stop,y); 
brv(x,stop+ 1 )=brv(x,stop); 
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brv(stop+ 1 ,stop+ 1 )=brv(stop,stop) j 

% 



Energy(e) 



bE(x,y)=oE(x,y)-dt*... 

((1/dx*... 

((oE(x+l ,y)+op(x+l ,y))*01ru(x+l ,y)/01r(x+l ,y) 
(oE(x,y)+op(x,y) )*01ru(x,y) /01r(x,y) )+... 
1/dy*... 

((oE(x,y+ 1 )+op(x,y+ 1 ))*01rv(x,y+ 1 )/01r(x,y+ 1 ) . 

(oE(x,y)+op(x,y) )*01rv(x,y) /01r(x,y) ))); 

bE(stop+l ,y)=bE(stop,y); 
bE(x,stop+ 1 )=bE(x,stop); 
bE(stop+ 1 ,stop+ 1 )=bE(stop,stop); 

% 



% return for next x 
end 

% return for next y 
end 

for x=2:stop+l 

for y=2:stop+l 



% 

% CORRECTOR 

% 

% Density(r) 

% 



nr (x,y)=(br(x,y)+01r(x,y))/2-... 

( dtdx/2*(bru(x,y)-bm(x-l ,y))+... 
dtdy/2*(brv(x,y)-brv(x,y- 1 ))); 



nr(l,y)=nr(2,y); 

nr(x,l)=nr(x,2); 

nr(2,2)=01r(2,2); 
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nr(l,l)=01r(2,2); 

Left wall boundary conditions 



% 



% Momentum(ru) 

% 

Px=gl*(bE(x,y)-0.5*((bru(x,y) A 2/br(x,y))+(bi-v(x,y) A 2/br(x,y)))); 

Pxm=g 1 *(bE(x- 1 ,y)-0. 5*((bru(x- 1 ,y) A 2/br(x- 1 ,y))+(brv(x- 1 ,y) A 2/br(x- 1 ,y)))); 

nru(x,y)=(bru(x,y)+01ru(x,y))/2-.. . 

((dtdx/2*... 

(( bru(x,y) *bru(x,y) /br(x,y) +Px )-... 

( bru(x- 1 ,y) *bru(x- 1 ,y)/br(x- 1 ,y)+Pxm)))+. . . 

(dtdy/2*... 

( bru(x,y) *brv(x,y) /br(x,y)-... 
bru(x,y- 1 ) *brv(x,y- 1 )/br(x,y- 1 ) ))) ; 



% 



nm(l,y)=nru(2,y); 

nru(x,l)=nru(x,2); 

nm(2,2)=01ru(2,2); 

nru(l,l)=01ru(2,2); 



% Momentum(rv) 

Py=gl*(bE(x,y)-0.5*((bru(x,y) A 2/br(x,y))+(biw(x,y) A 2/br(x,y)))); 

Pym=g 1 *(bE(x,y- 1 )-0. 5*((bru(x,y- 1 ) A 2/br(x,y- 1 ))+(brv(x,y- 1 ) A 2/br(x,y- 1 )))); 

nrv(x,y)=(brv(x,y)+01rv(x,y))/2-... 

((dtdx/2*... 

( bru(x,y) *bi"v(x,y) /br(x,y) -... 

bm(x- 1 ,y)*brv(x- 1 ,y)/br(x- 1 ,y) ))+. . . 

(dtdy/2*... 

( (brv(x,y)*brv(x,y)/br(x,y)+ Py )-... 

(brv(x,y- 1 ) *brv(x,y- 1 )/br(x,y- 1 )+Pym)))); 
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nrv( 1 ,y)=nrv(2,y); 
nrv(x,l)=nrv(x,2); 
nrv(2,2)=01rv(2,2); 
nrv(l,l)=01rv(2,2); 

% Left wall boundary conditions 



% Energy(E) 

% 



nE(x,y)=(bE(x,y)+oE(x,y))/2-... 

((dtdx/2*... 

((bm(x,y) /br(x,y) *(bE(x,y) +Px )) 

(bru(x- 1 ,y)/br(x- 1 ,y)*(bE(x- 1 ,y)+Pxm)))) +. . . 
(dtdy/2*... 

((brv(x,y) /bi-(x,y) *(bE(x,y)+Py )) 

(brv(x,y- 1 )/br(x,y- 1 ) *(bE(x,y- 1 )+Pym)))) ); 

% 

% 



nE(l,y)=nE(2,y); 

nE(x,l) = nE(x,2); 

nE(2,2)=oE(2,2); 

nE(l,l)=oE(2,2); 

%retum for next x 
end 

% 

%retum for next y 
end 

for y=l:stop+l 

for x=l:stop+l 

% 

% 

% 



Calculate pressure using the continuity equation for 
the ideal gas law and the new values for energy, 
momentum and density 
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np(x,y)=(g 1 ) *(n E(x,y)-(. 5 *((nru(x,y) A 2/. . . 
nr(x,y))+(nrv(x,y) A 2/nr(x,y))))); 



% — 
% 

% 

% 



Calculate velocity using the new values for momentum 
and density 



nu(x,y)=nru(x,y)/nr(x,y); 

nv(x,y)=nrv(x,y)/nr(x,y); 



% Check stability 

% Calculate the new speed of sound 

% 



nc=sqrt(g*np(x,y)/nr(x,y)); 

% 

% 

% Calculate the new time step dt 

% 



dtn=CFL*(l/(abs(nu(x,y))/dx+abs(nv(x,y))/dy+... 

nc*(l/dx A 2+l/dy A 2) A (0.5))); 
if dtn < dtmin 

dtmin=dtn; 

% end dtmin if statement 
end 



% 



% 



% 

%retum for next x 
end 

%retum for next y 
end 
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% 



% 

% Update time step and CFL(Courant number) 

dt=dtmin; 

fprintf('delta t= %1.5e s\tTotal time= %1.5e s\n',dt,tt) 



% 

% Output the values for pressure, density and velocity 

% using graphs 

% 

% Pressure 

for x=2:stop+l 

for y=2:stop+l 

press(x- 1 ,y- 1 )=np(x,y); 

end 

end 

colormap(cool) 
surf(press) 
hold on 

pcolor((press)-Pin) 

%contour(press,20) 
title('Absolute Pressure') 
xlabel('X grid points') 
ylabel('Y grid points') 
zlabel('Pascals’) 

pause(.Ol) 

%hold off 

% 

% Replace old values with new values and begin 

% next time step 

% 



op=np; 



Oh -nr; 
01ru=nru; 
01rv=nrv; 
oE=nE; 

% 

% 

end 

% 



Return for next t 
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